## Warning: package 'BiocGenerics' was built under R version 4.0.5
## Warning: package 'GenomeInfoDb' was built under R version 4.0.5
load("~/Documents/MiGASti/Databases/gene_matrix.RData")
metadata <- read.table("~/Documents/MiGASti/Databases/metadata2.txt")
#set rownames to Sample
row.names(metadata) <- metadata$Sample 
setwd("~/Documents/MiGASti/Databases")
exclude <- read.table("samples2remove2.txt")
exclude <- exclude$x
genes_tpm_filt = genes_tpm[, !colnames(genes_tpm) %in% exclude] 
#Excludes the samples from filters. 
#dim(genes_tpm_filt)
metadata_filt = metadata[ !(rownames(metadata) %in% exclude), ]
length(metadata_filt)
#check for expression of markers in dataset
markers = "~/Documents/MiGASti/Databases/Stimulations.xlsx"
markers = read_excel(markers, col_names = TRUE) 
markers = as.data.frame(markers)
gencode_30 = read.table("~/Documents/MiGASti/Databases/ens.geneid.gencode.v30")
colnames(gencode_30) = c("ensembl","symbol")
genes_tpm_filt = log2((genes_tpm_filt) + 1)
genes_tpm_filt <- as.data.frame(genes_tpm_filt)
setDT(genes_tpm_filt, keep.rownames = "ensembl")
genes_tpm_filt$ensembl
res_name = merge(genes_tpm_filt, gencode_30, by="ensembl")
rownames(res_name) = res_name$ensembl
marker_expression = merge(res_name, markers, by ="symbol")
dim (marker_expression)
head (marker_expression)
metadata_filt_ordered <- metadata_filt[colnames(genes_tpm_filt),]
metadata_filt_ordered = metadata_filt_ordered[-1,]
all(rownames(metadata_filt_ordered) == colnames (genes_tpm_filt)) #TRUEgit 
Stimulation <- metadata_filt_ordered$Stimulation

Expression of stimulation markers

rownames(marker_expression) = marker_expression$symbol
marker_expression$ensembl <- NULL
df_num = as.matrix(marker_expression[,2:484])
rownames(df_num) = marker_expression$symbol

#plot the genes based on samples 

metadata_filt$Stimulation <- as.factor(metadata_filt$Stimulation)
metadata_filt2 <- metadata_filt[order(metadata_filt$Stimulation, decreasing = T), ]
df <- t(df_num)
sorted <- df[rownames(metadata_filt2),]
sorted2 <- t(sorted)
sorted2 <- scale(sorted2)

ann <- data.frame(metadata_filt2$Stimulation)
colnames(ann) <- c('stimulation')
colours <- list('stimulation' = c('ununstim' = 'red2', 'unstim' = 'royalblue',  'TNFa'  = 'navy', 'R848'='green', 'LPS'='pink', 'IL4'='yellow', 'DEX'='orange', 'ATP'='brown', 'IFNy'='purple'))
colAnn <- HeatmapAnnotation(df = ann,
  which = 'col',
  col = colours,
  annotation_width = unit(c(1, 4), 'cm'),
  gap = unit(1, 'mm'))
#pdf(file = "~/Documents/MiGASti/Databases/Stim_markers.pdf", width = 100,
#height = 6)
colPalette <- colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(496)
Heatmap(sorted2,
        col = colPalette,
        name = "log2(TPM+1)", # Legend title
        cluster_rows = T,
        cluster_columns = F,
        #  row_names_gp = gpar(fontsize = 6), # Text size for row names)
        show_column_dend = F, 
        show_column_names = T,
        show_row_names = T,
        column_split = rep(c(metadata_filt2$Stimulation)),
        top_annotation = colAnn)
#dev.off()

#pdf(file = "~/Documents/MiGASti/Databases/Stim_markers.pdf", width = 100,
#height = 6)
colPalette <- colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(496)
Heatmap(sorted2,
        col = colPalette,
        name = "log2(TPM+1)", # Legend title
        cluster_rows = T,
        cluster_columns = F,
        #row_names_gp = gpar(fontsize = 4), # Text size for row names)
        show_column_dend = F, 
        show_column_names = T,
        show_row_names = T,
        column_split = rep(c(metadata_filt2$Stimulation)),
        top_annotation = colAnn)
#dev.off()

rownames(marker_expression) = marker_expression$symbol
marker_expression$ensembl <- NULL
df_num = as.matrix(marker_expression[,2:484])
rownames(df_num) = marker_expression$symbol

#plot the genes based on samples 

metadata_filt$Stimulation <- as.factor(metadata_filt$Stimulation)
metadata_filt2 <- metadata_filt[order(metadata_filt$Stimulation, decreasing = T), ]
df <- t(df_num)
sorted <- df[rownames(metadata_filt2),]
sorted <- scale(sorted)
sorted2 <- t(sorted)

Stimulation <- metadata_filt2$Stimulation

df <- rbind(sorted2, Stimulation)
library(dplyr)
df2 <- t(df)
df3 <- as.data.frame(df2)
setDT(df3)
df4 <- df3 %>% 
group_by(Stimulation) %>%
summarise_all("mean")

sorted_cor = as.matrix(df4[,2:27])
rownames(sorted_cor) <- c("ATP", "DEX", "IFNy", "IL4", "LPS", "R848", "TNFa", "unstim", "ununstim")
sorted_final <- t(sorted_cor)

cluster column = T, rows = F

# pdf(paste0(work_plots, "HM_markers_255s.pdf"), width = 10, height = 6)
colPalette <- colorRampPalette(rev(brewer.pal(5, "RdYlBu")))(15)
Heatmap(sorted_final,
        col = colPalette,
        name = "log2(TPM+1)", # Legend title
        cluster_rows = F,
        cluster_columns = T,
        #  row_names_gp = gpar(fontsize = 6), # Text size for row names)
        show_column_dend = T, 
        show_column_names = T,
        show_row_names = T)
#dev.off()

cluster row = T, columns = F

# pdf(paste0(work_plots, "HM_markers_255s.pdf"), width = 10, height = 6)
colPalette <- colorRampPalette(rev(brewer.pal(5, "RdYlBu")))(15)
Heatmap(sorted_final,
        col = colPalette,
        name = "log2(TPM+1)", # Legend title
        cluster_rows = T,
        cluster_columns = F,
        #  row_names_gp = gpar(fontsize = 6), # Text size for row names)
        show_column_dend = T, 
        show_column_names = T,
        show_row_names = T)
#dev.off()

cluster column and row = T

# pdf(paste0(work_plots, "HM_markers_255s.pdf"), width = 10, height = 6)
colPalette <- colorRampPalette(rev(brewer.pal(5, "RdYlBu")))(15)
Heatmap(sorted_final,
        col = colPalette,
        name = "log2(TPM+1)", # Legend title
        cluster_rows = T,
        cluster_columns = T,
        #  row_names_gp = gpar(fontsize = 6), # Text size for row names)
        show_column_dend = T, 
        show_column_names = T,
        show_row_names = T)
#dev.off()

Without ununstim

load("~/Documents/MiGASti/Databases/gene_matrix.RData")
metadata <- read.table("~/Documents/MiGASti/Databases/metadata.txt")
#set rownames to Sample
row.names(metadata) <- metadata$Sample 
setwd("~/Documents/MiGASti/Databases")
exclude <- read.table("samples2remove2.txt")
exclude <- exclude$x
genes_tpm_filt = genes_tpm[, !colnames(genes_tpm) %in% exclude] 
#Excludes the samples from filters. 
#dim(genes_tpm_filt)
metadata_filt = metadata[ !(rownames(metadata) %in% exclude), ]
length(metadata_filt)

[1] 38

#check for expression of markers in dataset
markers = "~/Documents/MiGASti/Databases/Signaling_pathways_2.xlsx"
markers = read_excel(markers, col_names = TRUE) 
markers = as.data.frame(markers)
gencode_30 = read.table("~/Documents/MiGASti/Databases/ens.geneid.gencode.v30")
colnames(gencode_30) = c("ensembl","symbol")
genes_tpm_filt = log2((genes_tpm_filt) + 1)
genes_tpm_filt <- as.data.frame(genes_tpm_filt)

#remove ununstim samples in metadata
metadata_cultured <- subset(metadata_filt, Stimulation == "unstim" | Stimulation == "ununstim")
#remove samples in genes counts datafile 
genes_tpm_cultured <- genes_tpm_filt[,metadata_cultured$Sample]

setDT(genes_tpm_cultured, keep.rownames = "ensembl")
res_name = merge(genes_tpm_cultured, gencode_30, by="ensembl")
rownames(res_name) = res_name$ensembl

marker_expression = merge(res_name, markers, by ="symbol")
rownames(marker_expression) = marker_expression$symbol
marker_expression$ensembl <- NULL
df_num = as.matrix(marker_expression[,2:172])
rownames(df_num) = marker_expression$symbol

#plot the genes based on samples 

metadata_cultured$Stimulation <- as.factor(metadata_cultured$Stimulation)
metadata_filt2 <- metadata_cultured[order(metadata_cultured$Stimulation, decreasing = T), ]
df <- t(df_num)
sorted <- df[rownames(metadata_filt2),]
sorted2 <- scale(sorted)
sorted2 <- t(sorted2)

ann <- data.frame(metadata_filt2$Stimulation)
colnames(ann) <- c('stimulation')
colours <- list('stimulation' = c('unstim' = 'royalblue',  'ununstim'  = 'navy'))
colAnn <- HeatmapAnnotation(df = ann,
  which = 'col',
  col = colours,
  annotation_width = unit(c(1, 4), 'cm'),
  gap = unit(1, 'mm'))


# pdf(paste0(work_plots, "HM_markers_255s.pdf"), width = 10, height = 6)
colPalette <- colorRampPalette(rev(brewer.pal(10, "RdYlBu")))(10)
Heatmap(sorted2,
        col = colPalette,
        name = "log2(TPM+1)", # Legend title
        cluster_rows = F,
        cluster_columns = F,
        #  row_names_gp = gpar(fontsize = 6), # Text size for row names)
        show_column_dend = F, 
        show_column_names = F,
        show_row_names = T,
        column_split = rep(c(metadata_filt2$Stimulation)),
        top_annotation = colAnn)